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ABSTRACT 

We describe a new set of self-consistent, equilibrium disk galaxy models that 
incorporate an exponential disk, a Hernquist model bulge, an NFW halo and a 
central supermassive black hole. The models are derived from explicit distribu- 
tion functions for each component and the large number of parameters permit 
detailed modeling of actual galaxies. We present techniques that use structural 
and kinematic data such as radial surface brightness profiles, rotation curves and 
bulge velocity dispersion profiles to find the best-fit models for the Milky Way 
and M31. Through N-body realizations of these models we explore their stability 
against the formation of bars. The models permit the study of a wide range of 
dynamical phenomenon with a high degree of realism. 

Subject headings: galaxies: individual (M31, Milky Way) — galaxies: structure 
— methods: N-body simulations — cosmology: dark matter 



1. Introduction 

The modeling of spiral galaxies to match photometric and dynamical measurements 
is a time-honoured endeavor. A complete and accurate characterization of the structural 
and kinematical properties of galaxies is necessary to understand the diversity of galactic 
dynamical behaviour and the origin of galaxies in the current cosmological paradigm. 
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This paper presents a new set of models for the phase-space distribution functions (DFs) 
of axisymmetric disk galaxies. The models consist of an exponential disk, a Hernquist model 
bulge (Hernquist 1990), an NFW halo (Navarro, Frenk k White 1996), and a central, super- 
massive black hole. They are defined by a large number of parameters which permit detailed 
modeling of real galaxies by fitting to observational constraints. Though the models repre- 
sent self-consistent equilibrium solutions to the coupled Poisson and collisionless Boltzmann 
(CB) equations, they are subject to both local and global nonaxisymmetric instabilities and 
are therefore suitable as initial conditions for numerical studies of galactic dynamics. 

The standard practice in modeling disk galaxies has been to assume simple functional 
forms for the space density and gravitational potential of the disk, bulge, and halo and 
then fit to a wide range of observational data. For the Milky Way galaxy, optimal structural 
parameters for these "mass models" are found primarily from surface brightness photometry, 
local stellar kinematics, the circular rotation curve, and observations of dynamical tracer 
populations such as globular clusters, the Magellanic clouds and the dwarf satellite galaxies. 
These observational constraints have been used to determine the best-fit model parameters 
based on ^-minimization techniques. (See, for example, Innanen (1973), Clutton-Brock, 
Innanen, k Papp (1977), Bahcall k Soneira (1980), Caldwell k Ostriker (1981), Kuijken k 
Gilmore (1991), Rohlfs k Kreitschmann (1988), Malhotra (1995); Kochanek (1996), Dehnen 
k Binney (1998a); Evans k Wilkinson (1999), Klypin, Zhao, k Sommerville (2002).) 

For external galaxies, the problem is somewhat simpler because of the advantage of 
an outside view and the smaller number of observational constraints. Optical and infrared 
de-projected surface brightness profiles of the stellar distribution are combined with rotation 
curves measured from HI gas kinematics to build mass models. The standard method is to 
perform a bulge-disk decomposition of the surface brightness profile usually into an i? 1//4 -law 
surface density profile for the bulge and a radial exponential profile for the disk. (See, for 
example, Kent (1985b), Simien, F. k de Vaucouleurs, G. (1986), and Binney k Merrifield 
(1998); see also Courteau, de Jong, k Broeils (1996) and references therein for a discussion 
of bulge-disk decomposition with the more general Sersic law for the bulge.) The halo model 
parameters are inferred by fitting multi-component mass models to the rotation curve data 
assuming values for the mass-to-light ratios of the disk and bulge stars. For nearby galaxies 
such as M31, kinematic data for the globular cluster, planetary nebulae and satellite systems 
are used to refine the models (Evans k Wilkinson 2000). 

There has also been parallel work in developing realistic N-body realizations of galaxy 
models for numerical experimentation in disk stability and galaxy interactions and mergers. 
To generate an N-body realization, it is necessary to model the full DFs for the multi- 
component system, in principle by solving the coupled Poisson and CB equations. Owing 
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to the complexity of these equations, various approximation schemes have been used. Early 
work on disk galaxy simulations set up initial conditions by estimating the velocity dispersion 
for a stable disk model with Toomre parameter Q > 1 and then assuming a locally Gaussian 
distribution for the velocities (see, for example, Sellwood (1985)). Barnes (1988) built 3- 
component models by adiabatically growing a disk potential within a live halo model and 
then adding an N-body disk with initial conditions generated as above. Hernquist (1993) 
presented a simpler approximate presciption for 3-component disk-bulge-halo models where 
the bulge and halo velocity distributions are assumed to be Maxwellian, truncated at the 
local escape speed and with dispersions estimated from the Jeans equations. These methods, 
while convenient, produce models that are slightly out of equilibrium and so require some 
relaxation time to damp out transients. Moreover, the models readjust to a state different 
from the one proposed (see, for example, Kazantzidis, Magorrian, & Moore (2004)). 

As numerical methods and resolution improve it has become necessary to develop more 
sophisticated techniques to generate initial conditions. Toward this end, Kuijken & Dubinski 
(1995, hereafter KD) presented a set of semi-analytic models for the DFs of disk galaxies 
consisting of an exponential disk, a centrally concentrated bulge, and an extended halo. 
The DFs are constructed from integrals of motion and represent equilibrium solutions to the 
Poisson and CB equations to a very good approximation and they have been used extensively 
to study different aspects of galaxy internal dynamics and interactions (e.g. Dubinski & 
Kuijken (1995), Dubinski (1998a), Dubinski et al. (1999), Garcia-Ruiz et al. (2002), O'Neill 
& Dubinski (2003), Widrow, Perrett, & Suyu (2003)). 

Fitting models to actual galaxies requires one to "observe" the model and compare 
with real observations. By providing the full DFs, the KD models enable one to add a 
level of realism to these pseudo-observations that is not possible with simple mass models. 
For example, a mass model rotation curve is constructed directly from the potential rather 
than from model line-of-sight velocities. With the KD models it is possible to construct 
the stellar rotation curve directly from stellar velocities and thereby incorporate asymmet- 
ric drift. Gravitational microlensing experiments provide another example where improved 
pseudo-observations are possible since the full DF allows for a self-consistent calculation of 
the predicted event rate distribution (Widrow, Perrett, & Suyu 2003). 

The bulge and halo of the KD models are characterized by a constant density core. Since 
their development it has become widely accepted from numerous cosmological simulations 
that dark halos have cuspy centers (Dubinski & Carlberg 1991; Navarro, Frenk & White 
1996). Detailed analysis of the central surface brightness and velocity dispersion profiles of 
early type galaxies and bulges suggest that these systems are also cuspy at their centres 
and furthermore contain supermassive black holes. We are therefore motivated to develop 
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a new set of axisymmetric disk galaxy models with cuspy halos and bulges that also allow 
for the self-consistent addition of a supermassive central object. These models have many 
purposes. They provide the gravitational potential and mass and velocity distributions from 
the sphere of influence of the central black hole in the inner few parsecs out to the virial radius 
of the model galaxy. Circular velocity curves, rotation curves incorporating asymmetric 
drift, velocity ellipsoids, line-of-sight velocity distributions (LOSVDs), and surface density 
(brightness) profiles can be generated easily. The large parameter space permits a wide range 
of models for comparison to real galaxies and the good quality of the initial equilibria make 
them ideal for studying subtle dynamical processes such as bar formation and disk warping 
and heating. 

The DFs are described in Section 2. Examples that match photometric and dynamical 
data for the Milky Way and M31 are presented in Section 3. This section also presents results 
from numerical experiments on the stability of these models as well as a brief discussion of 
an M31 model that incorportates a supermassive black hole. We conclude in Section 4 with 
a summary and a discussion of possible applications of the models. 

2. Distribution Functions 

The phase-space DFs for the disk, bulge, and halo of the KD models are chosen analytic 
functions of the integrals of motion. By Jeans theorem, any such DF yields a steady-state 
solution of the CBE in any potential that respects these integrals (Binney & Tremaine 1987). 
A self-consistent self-gravitating model is one in which the potential and space density also 
satisfy the Poisson equation. 

The KD models are, by design, axisymmetric with two known integrals of motion, the 
energy E and the angular momentum about the symmetry axis J z . KD choose the King 
model DF for the bulge which, in isolation, yields a system that is spherically symmetric and 
has a density profile characterized by a constant density core, an r~ 2 fall-off at intermediate 
radii, and a finite "tidal" radius where the density vanishes. For the halo, KD use the DF of 
a lowered Evans model which also exhibits a constant density core, power law intermediate 
region, and tidal radius. (We denote the tidal radius of the halo, which represents the outer 
edge of the system as a whole, as r t .) 

The configuration and velocity space distributions of the King and lowered Evans models 
are modified once they are incorporated into a multi-component model. Recall that for a 
system of collisionless particles in a static potential, any f(E) (i.e., any DF that is a function 
only of the energy) yields a steady-state solution of the CBE whose velocity distribution is 
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isotropic (Binney & Tremaine 1987). For an isolated self-gravitating system, p and \I> are 
necessarily spherically symmetric. (For a proof, see Perez & Aly (1996).) In the presence of 
an aspherical external potential, an f(E) will yield an aspherical mass distribution through 
the dependence of E on However, the velocities remain isotropic so long as / does not 
depend on angular momentum or some other integral of motion. For the multi-component 
models considered here and in Kuijken & Dubinski (1995), the disk potential causes a slight 
flattening of the bulge and halo. 

In this section we describe new models that have an NFW halo, a Hernquist bulge, and 
a central supermassive black hole. We begin by presenting DFs for the bulge and halo taken 
in isolation and then describe how these DFs are modified for composite models. 



2.1. The halo distribution function 

Navarro, Frenk & White (1996) found that the density profiles of dark matter halos in 
their cosmological simulations have a "universal" shape (the so-called NFW profile) of the 
form 

PNFw( r ) = — 2 (-0 

(r/a h )(l + r/a h ) 

where ah is the scale radius, ph = o\j^.-na\ is a characteristic density, and is a characteristic 
velocity dispersion. (Here and throughout we set Newton's constant G = 1.) In contrast 
with the profile of the Evans model halo, Pnfw has an inner r" 1 cusp and an extended r -3 
outer halo. The gravitation potential for this profile is 

2 log (1 + r/a h ) 

$NFW(?~) = ~<T h 1 • (2) 

r/a h 

Our strategy is to use a DF that, in the absence of a disk or bulge, yields a spherically 
symmetric NFW halo. We assume that the velocity dispersion is isotropic so that the DF 
depends only on the energy E. (Formally, E is the energy per unit mass.) For a given density 
profile, the DF can then be calculated through an Abel transform (Binney & Tremaine 1987) 
a procedure that has been carried out for the NFW profile by Zhao (1997), Widrow (2000), 
and Lokas & Mamon (2000). We write the DF as a function of the relative energy 8 = — E 



/nfw (£) = a h x a 2 JF NFW (S/al) (3) 
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For the dimensionless DF T we use the analytic fitting formula 

f , 3 e 3 / 2 (1 - er 5/2 f- log6 ) Q e p for < e < 1 
•W(e)= f 2 ^ [i e) ^^> e t0 " U<e<1 (4) 
j otherwise 

where P = P(£) is a fourth-order polynomial with P(l) = and ct = 2.71 (Widrow 2000). 

The NFW profile is infinite in extent and mass. For model building it is desireable to 
have a finite halo and so, following King (1966), we introduce an energy cutoff E h = —£ h = 
— tho\ and replace /nfw with 

/halo (£) = ^Ofc^halo (^A^) (5) 

where 

-r- ^ _ / -^nfw (e) - ^nfw (e/i) for e ft < e < 1 
^ (C) " \ other 



iU ; t/l . t " 1 (6) 
otherwise . 



e ft = corresponds to a full NFW profile while < e h < 1 yields a truncated profile. Exam- 
ples with ah — Oh — 1 and various values of are shown in Figure 1. Note that the truncation 
radius can be varied independent of the characteristic density and scale radius. When build- 
ing models a natural choice for the truncation radius is the cosmologically-motivated virial 
radius, r v - 1T (see below). After the inner properties of a model galaxy are set, a value of 
that gives r t m r vir can be easily found. In contrast, for the King and lowered Evans models 
a change in the truncation radius results in a significant change in the inner profile making 
model building more cumbersome and less intuitive. 

The DF as written above is symmetric under J z — > — J z and therefore generates a model 
with no net rotation. By spitting the DF into parts with positive and negative J z (J r ±) 
and recombining them with a suitable coefficient, one can generate a model with arbitrary 
amounts of rotation. Formally we write 

•Phalo = OL h J 7 + + (1 - OL h ) T- (7) 



where ah controls the amount of rotation (a^ = 1/2 implies no rotation). 
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2.2. The bulge distribution function 

Bulges are commonly modelled as the de Vaucouleurs r 1//4 law (de Vaucouleurs 1948) in 
projection and the Hernquist model is a simple density-potential pair which closely mimics 
this behaviour (Hernquist 1990). We therefore model the bulge using a Hernquist DF mod- 
ified by a energy cutoff E b to truncate the profile in a fashion similar to what was done for 
the halo. 

The standard Hernquist model has a density profile and potential given by 



Ph = 



Pb 



(r/a b ) (l + r/a b y 



(8) 



and 



-I 



1 + r/a b 

where a b , p b = al/2-Kal, and a b are the scale length, characteristic density, and characteristic 
velocity of the bulge. While the total mass is finite the density distribution is infinite in 
extent with p oc r~ 4 at large radii. 

We modify the Hernquist DF by incorporating an energy cutoff: 



/bulge (£) = v^afF^ (£ 1/2 M) (10) 

where 



T (n\ I T ^ ^ ~ T ^ (*) for g 6 < g < 1 

^ bulsc {q) ~ { otherwise, (U) 



q b = (—E/al) 1 ^ 2 , and JFH is the infinite-extent Hernquist model DF, 



MQ) = 27/^3 (^^5/2 ( 3 sin_1 9 + 9 C 1 - ^) C 1 - 2 (8^ 4 - 8g 2 - 3)) . (12) 
As with the halo, rotation is introduced through an additional parameter a b . 
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It is straightforward to use other models for the bulge. For example, the potential 
and DF for the density profile p oc r~ 3 / 2 (r + a)~ 5 ^ 2 are analytic and the associated surface 
density profile provides a somewhat better match to the de Vaucouleurs law than the surface 
density profile of the Hernquist model (Dehnen 1993). 

2.3. The disk distribution function 

The disk is assumed to be axisymmetric with space density paisk = Pdisk (R, z) and quasi- 
Maxwellian DF taken directly from the KD models which, in turn, were based on the 2-D 
model by Shu (1969) and extensions by Binney (1987). This DF is a function of E, J z , and 
E z , the latter being an approximate third integral of motion that corresponds to the energy 
in vertical oscillations. An implicit assumption in the formulation of this DF is that the 
velocity dispersions are small so that the epicyclic approximation is valid in the treatment 
of disk star orbits. This assumption limits the application of the models to cool disks. The 
DF can be tuned to yield a space density of a desired form. As in KD, we assume that the 
surface density profile of the disk is exponential in the radial direction with scale radius Rd 
and truncation radius i? ou t- The vertical structure is given approximately by sech 2 (z/zd) 
where Zd is the vertical scale height. In all, five parameters define the space density of the 
disk: R d , R ut, Zd, the mass M d , and a profile "shape" parameter 5R out which governs the 
sharpness of the truncation. 

The "stars" of the KD disks have nonzero velocity dispersions in the radial, azimuthal, 
and vertical directions. The dispersion in the radial direction, a R (R), is assumed to be 
exponential: cr R (R) = <7 R0 exp(—R/R a ). For simplicity, we set R„ = Rd in accord with 
observations by Bottema (1993). The dispersion in the azimuthal direction is related to a R 
through the epicycle equations (Binney & Tremaine 1987) while the dispersion in the vertical 
direction is set by the vertical potential gradient and the vertical scale height. 

In total, the DFs for the three components are characterized by fifteen parameters which 
we collect for convenience in Table 1. 



2.4. DFs for the composite model 

An isolated model halo is constructed by solving Poisson's equation self-consistently for 
the gravitational potential 



V 2 $ ha io = 47rp halo 



(13) 
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together with the integral relation between the DF (Eq.6), the density, and the relative 
potential ip = — $: 

Phaio = 2 5 / 2 tt d£^(£-i>)f halo (£) (14) 
Jo 

(Binney & Tremaine 1987). For a disk-bulge-halo model, Eq. 13 is replaced by 

V 2 $ = Anp (15) 

where 

p = J d 3 V (/disk + /bulge + /halo) (16) 

is the total density. 

The DFs in this expression depend on the energy (as well as other integrals of motion) 
which in turn depends on $. Since $ is the total gravitational potential, the density fields 
of the different components are inter-related in a complicated way. In fact the components 
of a disk-bulge-halo model constructed from Eq. 6 of KD and Eqs. 6 and 10 may bare little 
resemblance to the corresponding isolated components, a situation that is cumbersome for 
model building. 

To alleviate this problem, we proceed as follows. The DF for an isolated halo is nonzero 
over the energy range £ G {£h, cr%}- Qualitatively, we may expect the energy range of halo 
particles in the multi-component system to be extended to £ e {£h, o\ + o\ + where 
a\ = Md/Rd is the depth of the disk potential. These arguments suggest that we replace 
/halo {£) of Eq.5 by /h a io (£'h {£)) where £' h {£) is a function that maps the energy of a particle 
in the composite system {£) to the energy of a particle in the would-be isolated halo (£ f h )- 
In what follows, we assume 

£ = £k (i + °^ £ ') ■ ( 17 ) 

By construction, at low binding energies {£ — > 0), £' h ~ £ whereas £ — > a\ + of + o\ as 
£ 'h ~* °l- For the bul S e ' we re P lace /bulge {£) with / bulgc (£' b (£)) with £ b (£) = £' b + a 2 h + aj. 

These mappings are by no means unique and, apart from the conditions discussed above, 
are not motivated by any particular physical model. Indeed other mappings can be shown 
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to yield similar and equally acceptable models. 

In Figure 2 we show an example of the density profile and rotation curve for a typical 
disk-bulge-halo model. (The model parameters are given in Table 2.) Also shown is the halo 
profile that results assuming the same values of the at, o\ and Eh if the disk and bulge are 
not included. 



2.5. Comment on adiabatic compression 

It is important to stress that the "isolated halo" in Figure 2 is not meant to represent 
the progenitor for the halo in the final composite model. Indeed, the structure of dark halos 
in fully developed disk galaxies can be determined only through detailed modeling of galaxy 
formation. According to the standard structure formation paradigm at early times baryons 
and dark matter are well-mixed. Disks form when the (collisional) baryons lose energy (but 
not angular momentum) and settle to the bottom of the dark halo potential well. The halo 
in turn reponds to the change in the gravitational potential and therefore a pristine NFW 
halo will readjust to a new configuration. 

Simulations of galaxy formation which include gasdynamics and star formation are still 
relatively crude and computationally expensive. An alternative is to treat the effect of 
baryon infall on dark halos as an adiabatic process (Blumenthal et al. 1986; Flores et al. 
1993). If changes in the potential are slow compared with the orbital time of dark halo 
particles, then for each particle, there is a set of quantites known as adiabatic invariants 
which remain approximately constant even as the orbit changes. Knowledge of the adiabatic 
invariants allows one to determine the final DF from the initial DF without having to model 
the evolution of the system explicitly. 

Under a rather restrictive set of assumptions, the adiabatic theorem leads to the follow- 
ing simple relation: 

(Mbaryon(r) + M dm ) T = M ha i (ri)r; (18) 

(Young 1980; Blumenthal et al. 1986; Flores et al. 1993). In this formula, M baryon (r) is the 
baryon mass distribution as a function of radius (i.e., a spherically averaged representation 
of the disk and bulge), Md m is the final distribution of the dark halo, and Mh a io(^i) is the 
initial distribution of the baryon-dark matter proto-galaxy. Eq.18 allows one to go directly 
from the cosmologically-motivated NFW halo to the final dark halo in a fully-formed galaxy 
and has been used extensively in modeling disk galaxies (See, for example, (Mo, Mao, & 
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White 1998; Klypin, Zhao, & Sommerville 2002). 

In general, the dark halos predicted by Eq.18 are more centrally concentrated than the 
progenitors with a density profile that is cuspier in the center (e.g., r~ l cusp becomes an 
r~ L5 cusp as in Figure 7 of Klypin, Zhao, & Sommerville (2002)). However Eq.18 is based 
on a number of suspect assumptions. In addition to spherical symmetry and adiabaticity, 
one must assume that the dark halo particles are on circular orbits and therefore do not 
cross as their orbits contract. The particle orbits in simulated halos are, if anything, radially 
biased. (See Moore et al. (2004) who study the importance of the distribution of orbits for 
the evolution of dark halos.) 

An indication of the uncertainty in the adiabatic prescription can be seen in Figure 7 of 
Klypin, Zhao, & Sommerville (2002) where two halos from the same progenitor are shown. 
In one case, Eq.18 is used to determine the halo profile while in the other, baryons and dark 
matter are allowed to exchange angular momentum. The difference between the two halos 
is significant with the former being denser in the central regions by a factor of 30! 

How do our models fit in with the adiabatic compression paradigm? In short, we side- 
step the issue by focusing on models of fully developed disk galaxies. One might use Eq.18 
to work backwards from our final model and determine the progenitor (and thus compare 
with the NFW profile) but given the caveats associated with this equation we do not see 
this as a particularly fruitful endeavor. We do note that the difference between a model 
halo in isolation and one incorporated into a composite system is consistent with the general 
predictions of adiabatic compression theory — the profile exhibits a slightly cuspier inner 
region and is generally more concentrated — and so the naive use of the halo parameters as 
an indication of the general structure of the progenitor halo seems a reasonable first pass for 
making contact with cosmology. 

2.6. Models with central supermassive black holes 

The observation that most, if not all, disk galaxies harbor supermassive black holes near 
their centers has led to considerable interest in the structure of black hole-stellar systems 
as well as the effect a black hole might have on the central cusp of a dark matter halo. 
Tremaine et al. (1994) have derived DFs for a variety of models with central black holes 
whose density profiles have r~ v central cusps and r~ 4 outer parts. These so-called 77-models 
include the Hernquist profile (77 = 1) and are therefore directly applicable to the bulge DF 
in our system and easily generalized to the NFW halo model. The DFs of Tremaine et al. 
(1994) produce a density profile that is independent of the black hole mass: For a given 



77-model, the black hole alters the velocity distribution but not the space distribution of the 
stars and dark matter particles in its vicinity. We build black holes into our models under 
the same assumption. 

Consider first the bulge-black hole system. The effect of a black hole of mass M B h is to 
modify the relative potential 



^*(r) = ^(r) + — — (19) 
r 

where as in Tremaine et al. (1994) the superscript * denotes properties of models that include 
a black hole. While the Hernquist DF is non-zero over the energy range £ G {£ b , of}, with 
a black hole, the energy range is extended to £* G {£ b , oo}. 

In principle, disk-bulge-halo models with a central black hole can be constructed by 
replacing /bulge with a DF from Tremaine et al. (1994) that is modified by a suitable energy 
cut-off. The DFs for the disk and halo would be similarly replaced. However, the DFs in 
Tremaine et al. (1994) are not analytic and would have to be recalculated for each choice of 
A^bh- We therefore choose the following more efficient, albeit ad hoc scheme: The potential 
and density profile for a particular disk-bulge-halo model are calculated assuming no black 
hole. To incorporate a black hole, the potential is modified according to Eq.19 and a new DF 
is found which interpolates between /bulge of Eq. 10 and the DF of the appropriate ?7-model. 

At large binding energies (£* — > oo) the Hernquist-black hole DF has the asymptotic 
form (Tremaine et al. 1994) lim^.^oo /* = /^ where 



/c*\-V 2 



2 5 /% 3 M BH a fe W 
On the other hand, 

^2V^ 1 - £ ^T 5/2 PI) 

for £ — > of. From the asympototic form of the Hernquist potential, ipn — of / (1 — r/a b ), we 
have r ~ a b (1 + ipx/a^) suggesting 

£*-£ + M?L l \ (22) 



as a mapping from £ to £*. By construction, £*(£ b ) = £b- We take the DF to be 
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rr) = /^)(l + ^Py) 1 (23) 

which smoothly interpolates between /& and Eq.21. A similar procedure is carried out for 
the halo and disk. 



3. Models for the Milky Way and Andromeda Galaxies 

In this section we present models chosen to fit observational data for the Milky Way 
and Andromeda galaxies. The level of realism in our models is dictated to a large extent by 
the assumptions upon which they are based. At first glance, the assumption of axisymmetry 
seems rather restrictive since virtually all disk galaxies exhibit non-axisymmetric phenomena 
such as bars and spiral arms. However, our models are subject to non-axisymmetric instabil- 
ities. The program adopted here is to find the best-fit axisymmetric model and determine, 
through N-body simulations, if the model evolves to a state that more closely matches the 
actual galaxy. This approach has already been applied to the Milky Way by Sellwood (1985, 
1993), Fux (1997), and Valenzuela & Klypin (2003). 

Another key assumption is that the model galaxies are comprised of three components, 
a disk, bulge, and halo, whereas actual disk galaxies also have stellar halos and globular 
clusters that are typically spheroidal, more extended than the bulge, but (presumably) less 
extended than the dark halo. It would be straightforward to include such systems in our 
models but since they contain relatively little mass we do not do so here. The inclusion 
of globular cluster systems might be of interest for studying their evolution during galaxy 
mergers and interactions. 

Along similar lines, the structural parameters of the galactic disks (e.g., radial and 
vertical scale heights, velocity dispersion tensor) depend on colour and metalicity (Dehnen 
& Binney 1998b) whereas our models assume a single disk component. Stellar disks may, 
in fact, be more accurately represented as two component systems with a young thin disk 
and an old thick disk. Our models can be easily modified to include two or more disk-like 
components, an improvement which may prove relevant for detailed studies of bars and spiral 
structure. 

Our models assume DFs for the bulge and halo that are functions of the energy so 
that the velocity distribution in these components is necessarily isotropic. By contrast, the 
velocity distributions in simulated halos are typically biased toward radial orbits with a 
velocity anisotropy parameter (5 = 1 — Vq/v% ~ 0.6 (van den Bosch et al. 1999). Spherical 
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models with nonzero anisotropy parameters can be constructed from functions of E and 
J where J is the total angular momentum (Osipkov 1979; Merritt 1985; Baes & Dejonghe 
2002). However, J is not an integral of motion of a general axisymmetric. Moreover, since 
disks rotate, we can be fairly certain that halos do as well. Thus, realistic halo models require 
three or more integrals of motion. The velocity structure of dark matter particles will have 
an effect on the dynamical interaction between the halo and stellar components (e.g., the 
decay of the bar pattern speed through dynamical friction) and so it may be of interest to 
consider more general halo DFs. We leave the investigation of these subtle effects for future 
work. 

Perhaps the most severe assumption is that our models include only collisionless compo- 
nents whereas for some disk galaxies, a significant fraction of the "disk" mass is locked up in 
neutral hydrogen gas. Future versions of our models will include HI disks thereby expanding 
the applications to models with gas as well as stars. 

For each galaxy, we select a set of observational data which are compared with pseudo- 
observations of the model galaxy to yield a x 2 -statistic. Minimization of x 2 over the multi- 
dimensional parameter space yields the desired best-fit model. In addition, one can include 
non-observational constraints so as to select models with certain characteristics (e.g., specific 
value for the bulge-to-disk mass ratio, baryon fraction, or halo concentration parameter). 

In addition to the 15 parameters of Table 1 observer-dependent parameters such as the 
inclination angle for external galaxies or the galactocentric radius of the Sun for Milky Way 
must also be specified. In addition, if photometric data is used, then the mass-to-light ratios 
of the disk and bulge are also required. Depending on the type of observations, some of the 
parameters may be fixed during the minimization process. 

Following Widrow, Perrett, & Suyu (2003) we employ the downhill simplex algorithm 
(see, for example, Press et al. (1986)) for the minimization of \ 2 ■ A simplex is a geometrical 
figure in N dimensions where N is the number of parameters that defines the model. During 
execution of the algorithm, the simplex moves through parameter space seeking out the 
minimum of x 2 ■ As it does so, the simplex changes shape thus enabling it to move through 
complicated regions of parameter space. The downhill simplex method has a number of 
advantages over minimization procedures based on gradients of x 2 ( e -g- ; the method of 
steepest descent; see Press et al. (1986) and references therein). In particular, the method 
appears to be less susceptable to false minima though restarts are always executed to guard 
against this possibility. 
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3.1. The Milky Way 

Numerous authors have attempted to model the Milky Way with early important stud- 
ies by Bahcall & Soneira (1980) and Caldwell & Ostriker (1981). Bahcall & Soneira (1980) 
fit star count data together with local values for the scale heights and luminosity functions 
to parameterize global models of the disk and bulge. Caldwell & Ostriker (1981) consid- 
ered dynamical constraints (e.g., the Oort constants, rotation curve) in constructing three- 
component mass models. More recently, Dehnen & Binney (1998a) improved and updated 
these models by augmenting rotation curve data with dynamical constraints on the vertical 
structure of the Galaxy in the solar neighborhood. 



3.1.1. Observational constraints 

The models developed here for the Milky Way are assembled in the same spirit of early 
mass models but have the advantage that the end result is a fully realized DF for the stars 
and dark matter of the Galactic system. Seven observational data sets are used to constrain 
the Milky Way models. Five of these data sets, the inner and outer rotation curves, the 
Oort constants, the vertical force in the solar neighborhood, and the total mass at large 
radii are taken directly from Dehnen & Binney (1998a) and references therein. We also use 
measurements of the bulge dispersion at a projected distance of 200 pc from the galactic 
center (the peak of the dispersion profile) from the compilation of data by Tremaine et 
al. (2002). Finally, we incorporate estimates of the local velocity ellipsoid from Binney & 
Merrifield (1998). 

• Inner rotation curve 

Observations of HI emission provide a direct measure of the Galactic rotation curve. 
Inside the solar circle these observations are usually presented in terms of the so- 
called terminal velocity, t; term , the peak velocity along a given line-of-sight at galactic 
coordinates b = and |/| < 7r/2. Assuming the Galaxy is axisymmetric and the ISM 
rotates on circular orbits, the HI emission corresponding to t>t er m originates from the 
galactocentric radius R = R sml. Relative to the local standard of rest, we have 

vterm = v c (R) - v c (R ) sin I (24) 

where v c is the circular speed (see, for example, Binney & Merrifield (1998)). Following 
Dehnen & Binney (1998a) we use data from Malhotra (1995) restricted to the range 
sin / > 0.3 so as to avoid distortions from the bar. In particular, we use values of v teTm 
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from Figure 7 of her paper at four representative values of | sin/| averaging the results 
from the first and fourth quadrants. 

Outer rotation curve 

The radial velocity of an object at Galactic coordinates (/, b) relative to the local 
standard of rest, i>i sr is related to the circular rotation curve through the equation 

v isr = ( ~b-v c (R) - v c (R ) J cos 6 sin/ (25) 



, R 

1 /9 

where R = (d 2 cos 2 b + R$ — 2R d cos b sin /) and d is the distance to the object. 

In general, the data consists of a set of measurements (fi sr ,i, dj) which is to be compared 
with the model estimate W(R) and d(R) where R may be regared as a free parameter 
and W(R) = (Ro/R) v c (R) — v c (Ro) = vi sr / cos&sin/. For each data point, we adjust 
R so as to minimize 



Xi 



f W(R)-W i \ 2 ( d{R)-d l \ 2 

V AWi + { Ad t ) {2b) 



where W; L = v\ sri / cos b sin/. The x 2 f° r this data set is the average of the individual 
Xfs. 

In what follows we use data from Brand & Blitz (1993) with the same restrictions as 
in Dehnen & Binney (1998a) (i.e., / < 155° or / > 205°, d > 1 kpc, and W < 0). 

Vertical force above the disk 

Kuijken & Gilmore (1991) used K dwarf stars as tracers of the gravitational potential 
above the galactic plane thereby placing a constraint on the total mass surface density 
in the solar neighborhood. They found 



1^(1-1 kpc) 
2nG 



= 71 ± 6 M p C - 2 (27) 



independent of the relative contributions of the disk and halo. Only by including 
additional constraints on the local circular speed, galactocentric distance of the Sun, 
and Oort constants can one ferret out the separate contributions of the disk and halo 
to the local surface density. Doing so, Kuijken & Gilmore (1991) found 



E disk = 48 ± 9 M pc- 2 



(28) 
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which is in excellent agreement with estimates of known matter in the solar neighbor- 
hood. Since we include constraints on the Oort constants and rotation curve separately, 
we use Eq.27 as a constraint on our models and Eq.28 as a consistency check. 

• Oort constants 
The Oort constants 

and 

measure, respectively, the shear and vorticity in the Galactic disk. Following Dehnen 
& Binney (1998a), who review the published measurements, we adopt the constraints 

A = 14.5 ± 1.5kms _1 kpc _1 B = 12.5 ± 2 kms" 1 kpc' 1 . (31) 

• Local velocity ellipsoid 

The kinematics of stars in the solar neighborhood provides important constraints on 
the structure of the Milky Way. The observation that v 2 R ^ v 2 z already tells us that the 
disk DF cannot take the form f — f (E, L z ) and necessarily involves a third integral 
of motion (Binney & Tremaine 1987). Since the KD disks are built from three-integral 
DFs it is possible to model anisotropic velocity dispersion in galactic disks. There are 
two important caveats. The first, already mentioned in the introduction to this section, 
is that the models assume a single disk component so that the velocity dispersions in 
the radial, azimuthal, and vertical directions are single-valued functions of cylindrical 
radius. In point of fact, the shape of the velocity ellipsoid depends of colour (see, for 
example, Dehnen & Binney (1998b)). Furthermore, the velocity ellipsoid is rotated 
about the z-axis so that the principle axes do not coincide with the R and directions. 
The "vertex deviation" varies from 0° — 30° depending on B — V colour. 

These limitations can be overcome without too much difficulty (e.g., by including more 
than one disk component and by generalizing the disk DF to allow for vertex deviation 
of the velocity ellipsoid). We leave these improvements for future work and consider 
a single component disk. For constraints on the velocity ellipsoid, we use values from 
Table 10.4 of Binney & Merrifield (1998) which were derived from Edvardsson et al. 
(1993). Binney & Merrifield (1998) give values for the thin and thick disks and we use 



-18- 



a mass-weighted average assuming a 14:1 ratio between these two components (Dehnen 
& Binney 1998a) with 15% la errorbars. 

• Bulge dispersion 

Observations of the LOSVD in the direction of the bulge provide important constraints 
on the bulge parameters (and to a lesser extent, the parameters of the other compo- 
nents). Tremaine et al. (2002) have compiled measurements of the LOSVD between 
0.085 and 1300 pc. The dispersion profile shows a minimum of ~ 55 kms -1 at r ~ 5 pc 
and a maximum of 130 kms -1 at r ~ 200 pc (see, also, Kent (1992)). The rise of the 
dispersion profile inside r ~ 5 pc is presumably due to the central black hole while the 
detailed shape of the dispersion profile at larger radii may be affected by the barlike 
shape of the bulge. With this in mind, we average values of the line-of-sight dispersion 
near the peak to arrive at the single constraint cri os (R = 210 pc) = 136 ± 12kms _1 . 

• Mass at large radii 

The system of satellite galaxies that surround the Milky Way, the Magellanic Stream, 
and the high-velocity tail of the local stellar velocity distribution provide constraints 
on the large-scale mass distribution of the Galactic halo. Following Dehnen & Binney 
(1998a), who base their arguments on work by Kochanek (1996) and Lin, Jones, & 
Klemola (1995), we adopt 

M(r < 100 kpc) = (7 ±2.5) x 10 11 M (32) 
as a constraint on the mass distribution at large radii. 

3.1.2. Search strategy 

A x 2 statistic is calculated by comparing each of the seven data sets described above 
with pseudo-observations of the model. The pseudo-observations are designed to match 
closely the actual observations. For example, the LOSVD in the bulge region is found by 
calculating the velocity dispersion along a given line of sight of bulge "particles" chosen from 
the DF. In principle, one can add additional layers of realism to the pseudo-observations such 
as aperature smoothing for LOSVD measurements. 

The results are averaged in quadrature to yield a composite x 2 -statistic. In addition 
to the 15 parameters in Table 1 we must also specify the galactocentric radius of the Sun, 
Rq. We fix R out and 5R out to 30 kpc and 1 kpc respectively. Since the surface density of the 
disk falls exponentially, varying these values will not affect the model fit. In addition, the 
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rotation parameters of the halo and bulge are fixed so that neither component has net angular 
momentum. (In principle, incorporating more detailed observations of the bulge would allow 
us to fit the bulge rotation curve.) We pin the scale length of the radial dispersion profile to 
half that of the disk scale length (i.e., R a = so that the a\ and the surface density have 
the same exponential decay constant). We also fix = 0.2 which gives a tidal radius larger 
than 100 kpc. Finally, we run the simplex algorithm with R$ fixed to different values. In 
summary, each implementation of the simplex algorithm is run with nine free parameters. 

3.1.3. Results 

The parameter set for the best-fit Rq = 8 kpc model (MWa) is given in Table 2. Also 
given is the parameter set for a Milky Way model with a less massive disk that also has 
an acceptable x 2 (MWb). The density profiles and rotation curves for the two models are 
shown in Figures 3 and 4. We see that the disk dominates the rotation curve in model MWa 
for 3 kpc < R < 12 kpc whereas the disk never dominates the rotation curve of model MWb. 
We will return to this point below. 

A comparison of observations of tJ tcrm and W(R) with model predictions for MWa are 
shown in Figure 5 while a comparison of observed and predicted quantities for the other 
observables are given in Table 3. We see that both models provide an excellent fit to the 
observations. 

In Figure 6 we compare the line-of-sight velocity dispersion profile for model MWa with 
data compiled by Tremaine et al. (2002). Recall that to find the model, a single constraint at 
210 pc was used. We see that the agreement between model and observations is excellent for 
r > 100 pc. However, the model dispersion profile at smaller radii is too flat in comparison 
with the data. The slow rise in the dispersion profile appears to be a feature our models 
have in common with the 77-models of Tremaine et al. (1994) and may point to the necessity 
for a more complicated DF (e.g., one that depends on two or more integrals of motion). 

Figure 6 also includes the dispersion profile obtained for a model that includes a central 
black hole of mass 3.6 x 10 6 M Q . This value is roughly twice that found by Tremaine 
et al. (2002). The agreement between model and observations is poor between lpc and 
100 pc perhaps for the aforementioned reasons. However, the model profile does exhibit a 
minimum at r ~ 10 pc, a feature that is generic to galaxies that are believed to harbor 
central supermassive objects. 

How tightly do the data constrain the model? The answer in part is given in Figures 7 
and 8. For Figure 7, is fixed to different values while the remaining eight free parameters 
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are allowed to float. The result is x 2 as aa function of M d as shown in the upper panel. 
The derived disk-to-bulge mass ratio, again as a funtion of M d is shown in the lower panel. 
Evidently, the disk mass can vary by nearly a factor of two and still yield an acceptable fit 
to the data. 

A similar analysis for R is shown in Figure 8. Here, the minimum in x 2 (i? ) fairly 
flat especially toward larger values ot R . The preferred value of R ~ 8kpc agrees with 
well-known estimates. 

For both MWa and MWb, Rd/Ro — 0.36 which is approximately halfway between the 
two NFW models considered by Dehnen & Binney (1998a) (their models 2d and 4d). Table 
4 summarizes values for a number of other derived quantities for our models. Column 2 gives 
the bulge mass. The values for these two models are in excellent agreement with estimates 
of Mb based on COBE/DIRBE measurements (Dwek 1995) and gravitational microlensing 
observations (Bissantz 1997). 

Column 3 gives the Toomre Q parameter at 2.5Rd- Both models have Q > 1 indicated 
that they are stable against local perturbations (Binney & Tremaine 1987). Column 4 gives 
the local surface density of the disk. Our models evidently bracket the estimates from 
Kuijken & Gilmore (1991). 

Columns 5 and 6 give the total mass and tidal radius for the models. Recall that the tidal 
radius is controlled by the parameter e^. We have tuned this parameter so that r t is roughly 
equal to the virial radius i? v i r of the Galaxy as predected by the standard cosmological model 
of structure formation. In this scenario, -R v i r and the mass interior to this radius, M v - a , are 
related through the equation 

M vir = — A viT pR 3 yiT (33) 

where p is the mean density of the Universe and A vir parametrizes the average overdensity 
of the halo. (See Bullock et al. (2001) and references therein). For the currently favored 
ACDM cosmology (Q A = 0.3, Q m = 0.7) one finds A vir = 337. For both of our Milky Way 
models, r t ~ 240 kpc. 

Based on these values of r t , we can estimate the halo concentration parameter c v ; r ~ 
Tt/ah (Column 7). Model MWa compensates for the heavy disk by choosing a larger halo scale 
length and for this reason, its concentration parameter is smaller than that for MWb. Both 
models have concentration parameters larger than the mean value derived from cosmological 
simulations but within the la errorbars (Bullock et al. 2001). 
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3.1.4- Local dark matter density and the terrestrial dark matter detectors 

Dark matter detection experiments rely on halo models to develop search strategies and 
to analyse and interpret experimental data. Terrestrial dark matter detectors are sensitive 
to the local density and velocity-space distribution of dark matter particles. For these 
experiments, researchers have settled on a standard reference model, namely Maxwellian 
velocities with an rms speed of 270 km s -1 and a local density of p(Ro) = 0.0079 M pc~ 3 ~ 
0.3 (GeV/c 2 ) cm -3 . This reference model is useful for comparing the sensitivities of different 
experiments as well as making contact with predictions from theoretical particle physics. 

Our models allow one to study deviations from the standard model while ensuring that 
observational and dynamical constraints are satisfied. In Figure 9 we compare the local 
speed distribution for halo particles in models MWa and MWb with that of the standard 
reference model. The corresponding density is given in Column 8 of Table 4. We see that the 
heavy disk model (MWa) matches up quite well with the standard reference model while the 
light disk model has a local dark matter density that is a factor of two higher. As discussed 
above, MWa is unstable to bar formation suggesting a higher value for the local dark matter. 

Numerous authors have considered variations on the standard Galactic model such as 
bulk rotation of the halo, velocity space anisotropy, triaxiality, tidal streams, and small- 
scale dumpiness. Our models provide a starting point for further investigations along these 
lines. As described above, rotation may be added to the halo while velocity-space anisotropy 
and triaxiality require a non-trivial modification of the DF. In addition, our models are 
well-suited to numerical studies of the tidal disruption of subclumps in the dark halo. 

3.2. Models for M31 

We now take up the challenge of modeling M31 which provides a case study of an 
external galaxy. The surface brightness profile, the circular rotation curve, and the bulge 
velocity and dispersion profiles are combined to yield observation-driven models. 

Recently, Widrow, Perrett, & Suyu (2003) combined observations of the types described 
above to identify a suite of M31 models drawn from the original KD set. As one might expect, 
suitable models were found over a wide range in values for the disk mass-to-light ratio though 
models with particularly large values were found to be unstable to the formation of a strong 
bar and could therefore be ruled out. 

Our disk-bulge-halo DFs offer the possibility for an improved model of M31. The Hern- 
quist bulge is favored by observations while the NFW halo is favored by cosmological simu- 
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lations of structure formation. Morover, with the new models, the truncation radius of the 
NFW halo can be varied independent of the inner density profiles of the three components 
(Figure 1). Thus, one can tune the truncation radius to correspond to the virial radius as 
predicted by cosmology or alternatively, constrain the outer halo using observations of the 
M31 satellite system. 

3.2.1. Observations 

Following Kent (1989) and Widrow, Perrett, & Suyu (2003) we utilize measurements of 
Andromeda's surface brightness profile, rotation curve, and bulge velocity profiles. We use 
the global R-band surface brightness profile from Walterbos & Kennicutt (1987) which was 
obtained by averaging the light distribution over elliptical rings assuming an inclination of 
77°. 

The global surface brightness profile represents a small subset of the available photo- 
metric data. In principle a more sophisticated bulge-disk decomposition, where the thick- 
ness of the disk and inclination angle enter as free parameters, could be performed using 
two-dimensional surface brightness maps. However, apparent deviations of the galaxy from 
axisymmetry suggest that our models cannot adequately reflect this level of detail. For ex- 
ample, the position angle of the major axes of elliptical isophotes varies with galactocentric 
radius suggesting that the bulge is barlike and triaxial. Triaxiality can be introduced in 
a controlled way, for example by "adiabatically molding" the model into the desired shape 
(Holley-Bockelmann et al. 2001). Alternatively one can evolve, via N-body methods, a model 
galaxy that has a weak bar instability to see if some intermediate state is consistent with the 
observations (e.g., twisting isophotes). Such an exercise might further constrain the models 
and help break the mass-to-light degeneracy. 

For the rotation curve we combine measurements by Kent (1989) and Braun (1991) 
using the same kernel smoothing as in Widrow, Perrett, & Suyu (2003). The composite 
rotation curve extends from 2 — 25 kpc in galactocentric radii. Though the measurements 
in Braun (1991) extend to 30 kpc, we ignore data beyond 20 kpc since for this region of 
the galaxy, measurements were made along a single spiral arm on one side of the galaxy. 
Stellar rotation and velocity dispersion results from McElroy (1983) are used to constrain 
the dynamics of the inner 2 kpc of the galaxy. Widrow, Perrett, & Suyu (2003) attempted 
to fit the dispersion profiles along the major and minor axes as well as the bulge rotation 
profile between 300 and 2000 pc. Acceptable fits were found though in general, the models 
had a difficult time simultaneously reproducing a rising rotation curve and falling dispersion 
profile. One possible explanation is that the "bulge" data is contaminated by the disk (a 
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rapidly rotating, dynamically cold system). Here we use data between 300 and 1000 pc. 
The composite x 2 statistic is given by 

X 2 = ^(x!pb + xL + Xb) 1/2 ( 34 ) 

where Xspb> Xrc anc ^ xl f° r the surface brightness profile, rotation curve, and bulge bulk 
velocity and dispersion profiles. 



3.2.2. Search strategy 

Two of the seven parameters that describe the disk DF, the mass and radial scale length, 
are allowed to vary in the parameter-search algorithm. The vertical scale length, truncation 
radius, and truncation shape parameter are fixed at the values 0.3, 30, and 1 kpc respectively. 
Since our data do not depend on the dispersion of disk stars, the parameters am and R a 
are not required for the fitting algorithm though they are required for generating an N- 
body realization. The halo parameters and a h are allowed to vary while the truncation 
parameter is fixed to a value large enough so that the truncation radius is well outside the 
visible part of the galaxy. This parameter is adjusted after a suitable model for the visible 
part of the galaxy is found. Finally, the four parameters that describe the bulge, a&, a^, ft,, 
and ctb are allowed to vary. 



3.2.3. Results 

A wide range of models provide acceptable fits to the observations. Figure 10 is a 
contour plot of x 2 m the (M/L R ) d — (M/L R ) b plane ((M/L R ) d and (M/L R ) b are the ri- 
band mass-to-light ratios of the disk and bulge respectively). The general structure of a 
valley running approximately parallel to the (M/L R ) d -axis arises for two reasons. First, the 
bulge luminosity is constrained by the inner part of the surface brightness profile while its 
gravitational potential is constrained by the dispersion data. Hence, (M/L R ) b is relatively 
well-determined. By contrast, the primary constraint on the disk mass comes from the 
rotation curve. But since the disk and halo contributions to the rotation curve are similar, 
one can be played off the other and therefore the (M/L R ) d is poorly constrained. 

The degeneracy with respect to the disk mass-to-light ratio is generic to modeling of 
spiral galaxies and has been known for some time (see, for example, van Albada et al. (1985)). 
This degeneracy may be broken by fixing the disk mass-to-light ratio to a theoretically 
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preferred value, e.g., one derived from population synthesis models. Alternatively, one may 
constrain the disk mass-to-light ratio by requiring that the galaxy model be stable against 
the formation of a strong bar. As with the Milky Way, M31 may have a weak bar and so 
absolute stability is not a requirement (or even desirable). However, models with very heavy 
disks can clearly be ruled out as we demonstrate in the next section. A third possibility is 
to include an additional data set that probes more directly the disk mass distribution such 
as disk velocity dispersion measurements as in the study by Bottema (1987) of NGC5170. 

The mass-to-light ratios in Figure 10 must be corrected for foreground and external 
extinction if we are to make contact with theoretical predictions. Foreground extinction 
toward M31 is estimated to be 0.41 mag in B (de Vaucouleurs, de Vaucouleurs, & Corwin 
1976) corresponding to 0.23 mag in R assuming the standard interstellar extinction law 
(Binney & Merrifield 1998). Estimates of (R-band) internal extinction for the disk range 
from 0.6 (Monnet & Simien 1977) to 0.74 mag (Kent 1989) while formulae in Tully, et al. 
(1998) give 0.64 mag. Assuming 0.65 mag internal extinction for the disk and no internal 
extinction for the bulge (Kent 1989), (M/L R ) d and (M/L R ) b in Figure 10 should be scaled 
downward by 2.2 and 1.2 respectively to give x 2 in terms of intrinsic mass-to-light ratios. 

Population synthesis models provide an independent means of constraining the mass-to- 
light ratios in disk galaxies. Bell & de Jong (2001) provide a table of predicted mass-to-light 
ratios as a function of various color parameters. In particular, they find 

log 10 (M/L R ) = -0.820 + 0.851 (B - R) . (35) 

For the M31 disk, Walterbos & Kennicutt (1988) find B — R ~ 1.53 which, when corrected for 
extinction, corresponds to B — R ~ 1.18 implying (M/L R ) d ~ 1.5 (corrected) or (M/L R ) d ~ 
3.4 (uncorrected). Applying the same formula to the bulge region where B — R ~ 1.6 yields 
(M/L R ) b ~ 3.2 (corrected) or 3.9 (uncorrected). 

Our results may be compared with those from previous investigations. Kent (1989) 
modeled the disk and bulge of M31 using rotation curve and bulge dispersion measurements 
as well as photometric observations carried out in the r-bandpass of the Thuan & Gunn (1976) 
system. He found (uncorrected) values for the mass-to-light ratios of (M/L r ) disk = 10.5 and 
(M/L r ) bulgc = 5.0. Using the transformation r - R = 0.43 + 0.15 (B - V) from Kent (1985b) 
with a B — V color from Walterbos & Kennicutt (1988) Kent's values become (M/L R ) d ~ 6.2 
and (M/L R ) b ~ 3.0. Note that Kent (1989) assumes a constant density halo (i.e., halo with 
a very large core radius) which explains why his value for (M/L) d is higher than ours and 
others (Klypin, Zhao, & Sommerville 2002). 

Recent interest in bulge models has been driven in large part by an attempt to under- 
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stand their connection with galactic supermassive black holes. Magorrian et al. (1998) have 
constructed dynamical models for bulge-black hole systems in 36 spiral galaxies. For M31, 
they find a mass-to-light ratio in V-band of 4.83 ±0.1 or (M/L R ) b ~ 3.4 (uncorrected for 
extinction). 

To summarize, our analysis alone does little to constrain the disk mass-to-light ratio. 
As for the bulge, previous modeling as well as predictions from stellar population studies 
point to a mass-to-light ratio a factor of two higher than our preferred value. 

To better understand these results, we consider four models in more detail. Model M31a 
has disk and bulge mass-to-light ratios set equal to one another and to the value predicted by 
population synthesis models. This model lies close to the x 2 -valley in Figure 10. The same 
disk mass-to-light ratio is used in M31b but here (M/L R ) b is increased to 3.4 as estimated 
by Magorrian et al. (1998). The rotation curve and surface brightness profiles for these two 
models are shown in Figure 11 while the bulge line-of-sight dispersion and bulk velocity 
profiles are shown in Figure 12. The preference for lower values of (M/L R ) b is evident in 
the surface brightness and bulge line-of-sight dispersion profiles. With the higher value of 
(M/L R ) b , the predicted dispersion profile is too flat, i.e., does not fall fast enough with 
radius. The simplex algorithm strikes a balance between a bulge that is too dim and one 
that is too massive but ultimately, cannot fit the surface brightness and dispersion profiles 
as well as the low-(M/L R ) b model. Note that Kent (1989) has the same difficulty with the 
bulge dispersion profile while the analysis of Magorrian et al. (1998) is restricted to the 
innermost 0.15 arcmin (30 pc) of the bulge. 

M31c-M31a-M31d form a sequence of models from light to heavy disk mass. The rotation 
curves for models M31c and M31d are shown in Figure 13. 

3.3. Equilibrium and Stability 

We simulated N-body realizations of the four M31 models to test for both the quality 
of the derived equilibria and stability against bar formation. In all four models, the bulge 
dominates the innermost part of the rotation curve. In M31c, the disk contribution is 
subdominant throughout the system whereas in model M31d, the disk dominates the rotation 
curve between 5 and 25kpc. We anticipate that M31d will form a bar while M31c will be 
stable. For models M31a and M31b the disk contribution to the rotation curve is comparable 
to that of the other components at about one scale radius. 

For each model, we generate an N-body realization containing a total of 3.5M particles 
(1M disk particles, 500K bulge particles and 2M halo particles) and evolve the system using 
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a parallel N-body treecode (Dubinski 1998b). The extent to which the initial conditions 
represent a system in equilibrium and its susceptibility to the bar instability are explored 
by monitoring the surface density of the disk+bulge system, the disk velocity ellipsoid, the 
disk scale height, and the density profiles of the bulge and halo. We begin by discussing 
the results for model M31a which does not form a bar and represents a good test case. 
The central disk radial velocity dispersion has been set so that the Toomre Q parameter is 
Q ks 1.2 at it! = 2R d . This initial disk is rather cool and is unstable to spiral instabilities 
arising from the swing amplification of particle shot noise, an effect that leads to emerging 
spiral structure and disk heating. The disk is also heated through its interactions with halo 
particles but this effect is relatively minor with 2M halo particles and a halo particle mass 
is only rrih = 2.9 x 10 5 M . (A halo made of 10 6 M Q black holes will roughly double the disk 
scale height over a Hubble time (e.g. Lacey & Ostriker (1985)) and so our halo should have 
a much smaller effect.) 

We ran Model M31a for t = 4.8 Gyr using 10000 timesteps. The particle softening 
radius was set to s = 50 pc so the number of timesteps was more than adequate to follow 
orbits down to the softening radius in the galaxy core. The first impression is that the initial 
conditions are almost in perfect equilibrium with essentially no transient readjustment of 
the disk at start-up, i.e., no evidence of the imperfections in the initial conditions that are 
present in other methods (e.g. Hernquist (1993)). Figure 14 shows the disk+bulge surface 
density profile at t — 0, 2.4 and 4.8 Gyr along with the disk and bulge initial profiles. 
There is no change in the profile over the 4.8 Gyr integration and there is no sign of a bar. 
Figure 15 shows the evolution of the radially averaged disk velocity ellipsoid throughout the 
disk and reveals the effect of natural heating through the growth of spiral structure. The 
disk velocity ellipsoid at R = 15kpc (about 3 disk exponential scale lengths) changes from 
{(jRiCtfrCFz) = (22, 15, 18) km s _1 to (31,22, 19) km s _1 after 4.8 Gyr. Another indicator of 
disk evolution and heating is the vertical scale-height. Figure 16 shows the variance in disk 
particle height above the midplane as an estimator of disk scale-height at 3 times. There is 
some evolution with the scale-height growing by about 10% showing the good quality and 
relatively low shot noise effect in the models. The apparent flaring beyond about R = AR d 
may result partially from poor vertical force resolution due to discrete sampling of the disk. 
Since we choose a constant disk particle mass, the particle number density drops off quickly 
at large radii. 

Figure 17 shows the density profile of the bulge and halo over the course of the sim- 
ulation. The denser bulge develops a core radius of r c 300 kpc because of the use of 
force softening while the halo maintains its r _1 cusp. We see below in models with black 
holes that the cusp profile can be maintained if required using smaller softening lengths and 
smaller timesteps. We also examined the evolution of the halo shape profile as a further 
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test. Near the vicinity of the disk, the halo profile shape is flattened by the presence of the 
disk potential and so that flattening should remain unchanged if the system is in proper 
equilibrium. We compute the shape of the best fit ellipsoidal density contours as a function 
of radius at the start and end of the simulation using the normalized inertial tensor using the 
method described in Dubinski and Carlberg (1991) (Figure 18). This shows that the halo is 
flattened into an oblate spheroid with q ~ 0.8 near the disk and then becomes spherical at 
larger radii. With the exception of the central point where softening effects modify profiles, 
the shape profile remains unchanged throughout the simulation implying a good choice of 
equilibrium. In summary, the method produces an excellent equilibrium configuration of a 
spiral galaxy. The gradual heating of the disk can be attributed to the formation of transient 
spiral structure and to a lesser extent heating by halo particles. 

Simulations of the remaining models also show clean initial equilibria. Only model M31d 
develops a bar. The bar forms at 3 Gyr and persists until the end of the simulation at t — 4.8 
Gyr. At the end of the simulation the bar has a length of r& = 7.7 kpc with a pattern speed 
of 23 km s _1 kpc -1 and a co-rotation radius of Dl = 9.8 kpc. The ratio Dl : r& = 1.3 
makes this a "fast" bar in the standard nomenclature (e.g. (O'Neill & Dubinski 2003)) at 
least at this stage of the simulation. While M31 appears to have a triaxial bar-like bulge, 
at present, it does not appear to be a well-developed barred spiral. The absence of a bar in 
M31 suggests that the disk mass in Model D is too great and that one of the other models 
is more acceptable for M31. 

We also simulated the two Milky Way models and found little evolution in the velocity 
ellipsoid, vertical scale height, surface density, and space density again illustrating the excel- 
lent quality of the models as a means of setting up initial conditions for N-body simulations. 
Not surprisingly, model MWa (where the disk provides the dominant contribution to the 
rotation curve at intermediate radii) develops a bar while MWb appears to be stable. 

3.4. Incorporating a black hole 

In this section, we examine an axisymmetric model of M31 augmented by a central black 
hole. M31 is the nearest galaxy with a demonstrable central black hole but unfortunately 
exhibits a complex central structure in the form of a double nucleus (Lauer, et al. 1993; 
Bacon, Emsellem, Monnet, & Nieto 1994). Estimates of the black hole mass based on 
dynamical models and the observed bulge surface brightness and velocity dispersion profiles 
fall in the range 3 — 8.5 x 1O 7 M (Tremaine 1995; Kormendy & Bender 1999; Bacon et al. 
2001). 
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As an illustration of our models, we place a black hole of mass 3 x 10 7 M Q at the center 
of model M31a. Figures 19 and 20 show the line-of-sight velocity dispersion and surface 
density profiles. To generate the inner portion of the curves (dashed lines) in Figure 19 we 
sample 10M particles from the bulge DF within a thin tube 20 pc in radius centered on the 
galaxy and aligned with the symmetry axis. The dispersion profile exhibits the characteristic 
intermediate minimum and maximum seen in the data (Tremaine 1995). 

We also simulate a galaxy model including a blackhole to test the quality of the equi- 
librium focusing on the nuclear region close to the blackhole. The blackhole moves freely 
as single particle within the simulation. For this test, we put most of the particles in the 
bulge, iV = 10M so the ratio of the blackhole mass to the bulge particle mass is ~ 10 4 . The 
enclosed mass of bulge stars equals the blackhole mass at a radius r = 34 pc. there are only 
The halo and disk are represented by 2M and 1M particles respectively but their particles 
are unlikely to interact with the blackhole in this test. We set the particle and blackhole 
softening length to s = 1 pc and use a single timestep of At = 100 yrs. Even with a total 
of N = 10M particles in the bulge, there are only 13 bulge particles initially within the 
softening radius. Figure 20 also shows the evolution of the model focusing on the central 
region surrounding the black hole. The simulation is run for 3000 timesteps corresponding 
to approximately 18 orbital times at a radius of 1 pc. We see that there is some evolution 
of the system — towards the center, the velocity dispersion fluctuates slightly within r ~ 10 
pc and the surface density increases slightly. At r = 1 pc where there are only a dozen or 
so particles, we expect larger initial variations that may explain the settling to a slightly 
higher density. Nevertheless, the profiles remain close to the initial state suggesting that 
the equilibrium is reasonable. These effects may be also result indirectly from using force 
softening or a reflection of the approximations that have gone into the DF. Clearly, high 
numerical resolution is needed to treat the dynamics of stars around a central blackhole and 
the short orbital times make these technically difficult and costly. 

Models such as this one may be used as initial conditions for numerical experiments of 
black hole dynamics during galaxy merger events such as the simulations by Milosavljevic, 
& Merritt (2003) and will be explored in future investigations. 

4. Summary and Discussion 

Our goals in this paper have been threefold. First we have presented a new set of model 
DFs for multi-component disk galaxies. Second, we have identified particular models that 
fit observational data for the Milky Way and M31. Finally, we have explored the stability 
of the models using numerical simulations. 
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Our DFs respresent self-consistent axisymmetric equilibrium solutions to the Poisson 
and CB equations. The disk and bulge DFs are motivated by observations while the halo 
DF is motivated by results from cosmological N-body simulations. The models permit the 
inclusion of a central supermassive black hole. 

Historically, two approaches have been used to construct DFs for equilibrium systems 
(Binney & Tremaine 1987). The first approach is to propose a DF built from chosen functions 
of the integrals of motion and then calculate the density by solving the Poisson equation. 
The second is to calculate the DF from the desired density profile via an integral equation 
known as the Abel transform. Our DFs are the product of a hybrid scheme. The DF for the 
isolated NFW halo was derived from the Abel transform, modified by the energy mapping 
as described in the text, and then used for the halo in the composite system. 

Our models are defined by 15 free parameters which may be tuned to fit a wide range of 
observational data. Typically most of the parameters are poorly constrained. A galaxy's sur- 
face brightness profile is, in general, adequate to fix the disk scale length and, in combination 
with velocity dispersion and rotation curve measurements, provides a constraint on the disk 
mass and halo and bulge structural parameters. However degeneracies remain primarily due 
to uncertainties in the disk and bulge mass-to-light ratios. Theoretical considerations such 
as results from population synthesis studies or numerical experiments of the bar instability 
can help narrow down the field of acceptable models. Additional data for external galaxies 
such as velocity dispersion measurements in the disk may help break the M/L degeneracy. 
Observations of edge-on spiral galaxies can be used to constrain the vertical scale height and 
disk trucation parameters (Kregel, van der Kruit, & de Grijs 2002). 

We have constructed sequences of models for the Milky Way and M31 which provide 
excellent fits to available data. The models serve to illustrate the general proceduce for 
searching a large parameter space to find acceptable models for particular galaxies. 

A primary purpose of our models is to provide initial conditions for N-body experi- 
ments and well-developed techniques allow one to sample the DFs with arbitrary numbers 
of particles. Through a series of numerical experiments, we have explored the quality of the 
models as initial conditions and their stability to the formation of bars. Models in which the 
disk contribution to the rotation curve is always subdominant to that of the bulge and halo 
tend to be stable against bar formation and therefore provide the best laboratory to study 
the quality of the model DFs. Our analysis of one such model indicates that the DFs do 
indeed provide excellent initial conditions. The radial profiles of the space density, surface 
density, velocity dispersion tensor, and vertical scale height all remain relatively constant 
over 4.8 Gyr. Areas where we do see some evolution may be understood from simple con- 
siderations. For example, the bulge evolves from Hernquist cusp to core inside the softening 
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length because the velocity distribution is calculated assuming an unsoftened force-law. Spi- 
ral density waves, generated by swing-amplification of shot noise cause the disk to thicken, 
but only by about 10%. 

Models in which the disk contribution to the rotation curve is dominant over some range 
in radius, develop strong bars. We confirm this well-known result for both Milky Way and 
M31 models. 

Our models allow for the addition of a central supermassive black hole. Our prescription 
for doing this maintains the density distribution of the black hole-less models while modifying 
the velocities to establish a new dynamical equilibria near the black hole. The velocity 
dispersion profiles of our black hole models have the right characteristics to match the data, 
namely high dispersion near the black hole falling to a minimum at tens of parsecs, rising at 
a radius of a hundred or so parsecs, and then falling again as one moves further out in radius. 
Thus, we may be able to model galaxies from the sphere of influence of the black hole out 
to the virial radius, an impressive dynamic range of five or more orders of magnitude. 

There are numerous open avenues for future work. The model fitting algorithm can ac- 
comodate any number and type of observations and one may add to the pseudo-observations 
layers of realism such as statistical fluctuations and the effects of seeing. Our M31 study used 
a small subset of the available photometric and kinematic data. Two-dimensional surface 
brightness and dispersion maps reveal the presence of the bar and may be compared with 
corresponding maps from the "evolved" models providing a further constraint on the model. 

A classic problem in galactic dynamics is the disruption of satellite systems by the tidal 
field of the parent galaxy. The discovery of an arc of stars associated with the Sagittarius 
dwarf galaxy (see Majewski et al. (2003) and references therein) has intensified interest in 
performing simulations of the tidal disruption of this system (Ibata et al. 2001; Helmi 2004) 
(see also Geehan, Fardal, Babul, & Guhathakurta (2005) for a study of the Andromeda 
stream in a disk-bulge-halo model of M31). Our models offer the possibility of performing 
such simulations realistic and fully-self-consistent models for both parent and satellite sys- 
tems. Since a primary goal of previous investigations has been to constrain the shape of the 
Galactic halo, our models will have to be extended to include triaxial systems if they are 
to be useful in this endeavor. The extension to triaxial systems might be accomplished by 
adiabatically molding the models as in Holley-Bockelmann et al. (2001). 

Our models might be extended in other important ways such as the inclusion of a 
globular cluster system or thick disk as discussed in the text. Perhaps the most significant 
and challenging improvement would be to add a gas component and star formation since 
this would enable the study of spiral structure and bar formation with a far greater level of 
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realism. 
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Fig. 1. — Density profiles (top) and rotation curves (bottom) for halo models with ah = 
ah — 1 and values of Eh corresponding to tidal radii r t = 5, 10, 40, and 80. Also shown is the 
full NFW model (solid curve). 
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Fig. 2. — Density profile (top) and rotation curve (bottom) for the disk-bulge-halo model 
described in the text. Shown are the contributions from the disk (dotted curve), bulge 
(dashed curve), halo (long-dashed curve) as well as the total density profile and rotation 
curve (solid curve). Also shown is the halo model that results if the same halo parameters 
are used and the disk and bulge are not included (dot-dashed curve in the top panel). 
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Fig. 3. — Density profile and rotation curve for model MWa. Line types are the same as in 
Figure 2. 




Fig. 4. — Density profile and rotation curve for model MWb. 
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Fig. 5. — Results from the best-fit Milky Way model compared with observations for W 
(top) and fterm (bottom) as defined in the text. 
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Fig. 6. — Line of sight dispersion profile for the bulge region. Data are from a compilation 
by Tremaine et al. (2002) of published measurements. 
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Fig. 7. — x 2 -statistic as a function of M d (top). The best-fit model is found with M d fit. For 
each case, we show the bulge mass M b as a function of M d . 
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Fig. 9. — Speed distribution of halo particles at the position of the solar system for model 
M31a (dotted curve), M31b (dashed curve), and the standard model used by experimentalists 
(solid curve). 
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Fig. 11. — Predictions vs. observations for models M31a (left) and M31b (right). Shown are 
the rotation curve (top), R-band surface brightness profile (middle) and residuals (predicted 
- observed) for surface brightness profile (bottom). Linetypes in the rotation curve plots are 
the same as in previous figures. 
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Fig. 12. — Bulge dispersion and bulk rotation profiles for models M31a (solid curve), M31b 
(dotted curve) and observations (data points). 
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Fig. 13. — Rotation curves for models M31c (left) and M31d (right). 
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Fig. 14. — Evolution of the surface density profile over 4.8 Gyr. The three solid curves give 
the total surface density at t — 0, 2.4, and 4.8 Gyr. The dotted curves show the separate 
contributions from the bulge and disk. 



-51 - 




10 20 30 

R (kpc) 



Fig. 15. — Evolution of the disk velocity ellipsoid over 10 Gyrs. The profiles of the velocity 
dispersions are shown at t — 0.0, 2.4 and 4.8 Gyr showing a slow increase over the disk. 
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Fig. 16. — Disk scale height evolution at t — (solid curve) and t = 2.4 Gyr and t = 4.8 
Gyr (dotted curves) 



-53- 




-2-10123 



log R/kpc 



Fig. 17. — Bulge and halo radial density profiles over 4.8 Gyr. The vertical dotted line 
corresponds to the softening length. A constant density core develops in the bulge on the 
scale of the softening length of s = 50 pc. 
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Fig. 18. — Halo axis ratio profile at the start (solid) and end (dashed) of the simulation over 
4.8 Gyr. With the exception of the central point the axis ratio of the halo remains essentially 
constant throughout the simulation. 
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Fig. 19. — Velocity dispersion (top) and surface brightness (bottom) profiles in the bulge. 
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Fig. 20. — Evolution of the velocity dispersion and surface brightness profiles around the 
central black hole. Solid curve - t — 0; dotted curve - t ~ 10 5 yrs short-dashed curve - 
t ~ 2 x 10 5 yrs long-dashed curve - t~3x 10 5 yrs. The orbital time for the black hole at 
lpc is 1.7 x 10 4 yrs. 
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Table 1. summary of model parameters 



parameter description 





halo tidal radius parameter 


&h 


halo characteristic velocity 


ah 


halo scale length 


ah 


halo bulk rotation 


M d 


disk mass 


R d 


disk scale length 


Rout 


disk truncation radius 


SRout 


sharpness of truncation 


h d 


disk scale height 


oro 


radial velocity dispersion at galaxy center 


Ra 


scale length for radial dispersion 




bulge tidal radius parameter 


Ob 


characteristic bulge velocity 


a b 


bulge scale length 


otb 


bulge rotation 



Table 2. parameters for models discussed in the text 



Model 


th 


Oh 


ah 


M d 


Rd 


hd 


£6 


Ob 


ab 




(M/L R ) d 


(M/L R ) h 


Sample 


0.079 


1 


1 


0.1 


0.3 


0.02 


.1 


1.15 


0.15 








MWa 


0.17 


2.496 


12.96 


19.66 


2.806 


0.409 


0.787 


4.444 


0.788 


1.211 






MWb 


0.11 


3.447 


8.818 


14.47 


2.817 


0.439 


0.791 


4.357 


0.884 


1.244 






M31a 


0.25 


3.371 


12.94 


33.40 


5.577 


0.3 


0.929 


4.607 


1.826 


0.763 


3.4 


1.9 


M31b 


0.25 


3.224 


14.03 


35.08 


5.401 


0.3 


0.925 


4.811 


1.857 


0.751 


3.4 


3.4 


M31d 


0.36 


3.243 


17.46 


50.10 


5.566 


0.3 


0.924 


4.685 


1.802 


0.732 


5.0 


2.5 



Note. - 



We assume G = 1. Units for the Milky Way and M31 models are kpc, 100 kms" 1 and 2.33 x 10 9 M< 
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Table 3. Comparison of observations and model predictions 

























x 2 


A 


B 


< 


v% 


K z /2irG 


(7los(210pc) 


M(100kpc) 


observed 




14.5 


-12.5 


36 


25 


20 


71 


136 


70 


MWa 


0.75 


13.6 


-12.5 


31.7 


26.4 


19.9 


71.7 


134 


67 


MWb 


0.91 


13.5 


-12.9 


32.3 


27.5 


20.0 


68.6 


134 


72 



Note. — Units are kms 1 for the velocity dispersions, 10 10 M Q for Mb and M(100kpc) and 
M pc- 2 for K z /2nG. 



Table 4. derived quantities for Milky Way models 



Model 


M h 


Q 


Kz.yis/2-KG 


M t 


n 


Ct 


p(Ro) 


MWa 


1.1 


1.3 


53 


74 


239 


19 


0.0079 


MWb 


1.2 


2.2 


39 


77 


236 


27 


0.015 



Note. — Units are kpc for r t and M Q pc 3 for p{Rq). Other- 
wise the units are the same as in Table 2. 



